Aeolian transport layer 
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We investigate the airborne transport of particles on a granular surface by the saltation mechanism 
through numerical simulation of particle motion coupled with turbulent flow. We determine the 
saturated flux q 3 and show that its behavior is consistent with a classical empirical relation obtained 
from wind tunnel measurements. Our results also allow to propose a new relation valid for small 
fluxes, namely, q s = a(u* — ut) a , where u, and ut are the shear and threshold velocities of the wind, 
respectively, and the scaling exponent is a ~ 2. We obtain an expression for the velocity profile of 
the wind distorted by the particle motion and present a dynamical scaling relation. We also find 
a novel expression for the dependence of the height of the saltation layer as function of the wind 
velocity. 

PACS numbers: 45.70.Mg, 47.55.Kf, 47.27.-i, 83.80.Hj 
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The transport of sand by wind is among others re- 
sponsible for sand encroachment, dune motion and the 
formation of coastal and desert landscapes. The domi- 
nating transport mechanism is saltation as first described 
by Bagnold 1] which consists of grains being ejected up- 
wards, accelerated by the wind and finally impacting onto 
the ground producing a splash of new ejected particles. 
Reviews are given in Refs. (EH- Quantitatively this pro- 
cess is however far from being understood. 

Due to Newton's second law the wind loses more mo- 
mentum with increasing number of airborne particles un- 
til a saturation is reached. The maximum number of 
grains a wind of given strength can carry through a unit 
area per unit time defines the saturated flux of sand q s . 
This quantity has been measured by many authors in 
wind tunnel experiments and on the field, and numerous 
empirical expressions for its dependence on the strength 
of the wind have been proposed 0, IE IE ■ I n previ- 
ous studies theoretical forms have also been derived us- 
ing approximations for the drag in turbulent flow [E ■ 
All these relations are expressed as polynomials in the 
wind shear velocity which are of third order, under 
the assumption that the grain hopping length scales with 
and otherwise can be more complex 
@. The velocity profile in a particle laden layer has also 
been the object of measurement [Tl], [l2| and modelliza- 
tion [l3| . Surprisingly however very few measurements of 
the height of the saltation layers as function of it* have 
been reported [T^ | and no systematic data close to the 
threshold are available. The complete analytical treat- 
ment of this problem remains out of reach not only be- 
cause of the turbulent character of the wind, but also 
due to the underlying moving boundary conditions in the 
equations of motion. Despite much research in the past 
there remain many uncertainties about the trajecto- 
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FIG. 1: Schematic representation of the setup showing the 
mobile wall at the top, the velocity field at different posi- 
tions in the ^-direction and the trajectory of a particle stream 
(dashed line). At the collision between this stream and the 
static wall at the bottom, we consider that the particles re- 
bound to the air flow at a ejection angle 6, with only a fraction 
r of their original kinetic impact energy. 



ries of the particles and their feedback with the velocity 
field of the wind. It is this challenge which motivates the 
present work. 

We will present the first numerical study of saltation 
which solves the turbulent wind field and its feedback 
with the dragged particles. As compared to real data, our 
values have no experimental fluctuations neither in the 
wind field nor in the particle sizes. As a consequence, we 
can determine all quantities with higher precision than 
ever before, and therefore with a better resolution close to 
the critical velocity at which the saltation process starts. 

In order to get quantitative understanding of the layer 
of airborne particle transport above a granular surface, 
we simulate the situation inside a two-dimensional chan- 
nel with a mobile top wall as schematically shown in 
Fig. 1. We impose a pressure gradient between the left 
and the right side. Gravity points down, i.e., in negative 
^/-direction. The ^-dependence of the pressure drop is ad- 
justed in such a way as to insure a logarithmic velocity 
profile along the entire channel in the case without parti- 
cles, as it is expected in fully developed turbulence [Tt3 |. 
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More precisely, this profile follows the classical form 

u x (y) = (u*/K)ln(y/yo) , (1) 

where u x is the component of the wind velocity in the 
x-direction, it* is the shear velocity, k = 0.4 is the von 
Karman constant and yo is the roughness length which 
is typically between 10~ 4 and 10~ 2 to. The upper wall of 
the channel is moved with a velocity equal to the velocity 
of the wind at that height in order to insure a non-slip 
boundary condition. 

The fluid mechanics inside the channel is based on the 
assumptions that we have an incompressible and Newto- 
nian fluid flowing under steady-state and homogeneous 
turbulent conditions. The fluid is air with viscosity p = 
1.7894 x I0~ 5 kg m^s^ 1 and density p = 1.225 kg m~ 3 . 
The Reynolds-averaged Navier-Stokes equations with the 
standard k — e model are used to describe turbulence. 
The numerical solution for the velocity and pressure fields 
is obtained through discretization by means of the con- 
trol volume finite-difference technique 0, 0] . The inte- 
gral version of the governing equations is considered at 
each cell of the numerical grid to generate a set of non- 
linear algebraic equations which are pseudo-linearized 
and solved. The criteria for convergence used in the sim- 
ulations are defined in terms of residuals, i.e., a measure 
of the degree to which the conservation equations are sat- 
isfied throughout the flow field. In our simulations con- 
vergence is achieved only when each of the normalized 
residuals falls below 10~ 6 . 

After having produced a steady-state turbulent flow, 
we proceed with the simulation of the particle transport 
along the channel. Assuming that drag and gravity are 
the only relevant forces acting on the particles, their tra- 
jectory can be obtained by integrating the following equa- 
tion of motion: 

^ = F D (u-u p )+g(p p -p)/p p , (2) 

where u p is the particle velocity, g is gravity and p p — 
2650 kg to~ 3 is a typical value for the density of sand 
particles. The term Fd(u — u p ) represents the drag force 
per unit particle mass where 
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dp = 2.5 x 10~ 4 rn is a typical particle diameter, Re = 
pd p \u p — u| / p is the particle Reynolds number, and the 
drag coefficient Co is taken from empirical relations 
|19| . Each particle in our calculation represents in fact a 
stream of real grains. It is necessary to take into account 
the feedback on the local fluid velocity due to the mo- 
mentum transfer to and from the particles. Specifically, 
this coupling effect is considered here by alternately solv- 
ing the discrete and continuous phase equations until the 
solutions in both phases agree. The momentum transfer 
from one phase to another is computed by adding the 
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FIG. 2: Typical trajectories of particles computed for u* < u t 
(full line) and u* > u t (dashed line). 



momentum chan ge o f every particle as it passes through 
a control volume llsTl. 



F D (u- u p )rh p At 

particles 
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where m p is the mass flow rate of the particles and At 
the time step. The exchange term Eq. Q appears as a 
sink in the continuous phase momentum balance. 

In Fig. 1 we see the trajectory of one particle stream 
and the velocity vectors along the ?/-direction. Each time 
a particle hits the ground it loses a fraction r of its energy 
and a new stream of particles is ejected at that position 
with an angle 9. The parameters r — 0.84 and = 36° 
are chosen from experimental measurements [53, 0] • We 
also studied other values for r and 9 and even considered 
a continuous distribution of ejection angles. As expected, 
the choice of unrealistic values produces unphysical re- 
sults. More details will be given in Ref. 22]. 

If it* is below a threshold value Ut the energy loss at 
each impact prevails over the energy gain during the ac- 
celeration through drag and particle transport comes to 
a halt as illustrated in Fig. 2. Only for it* > it t steady 
sand motion is achieved. The resulting flux is given by 
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where m p and (u p ) are the mass and the average velocity 
of the particles, respectively, and n p is the number of par- 
ticle streams released. The first added particle streams 
are strongly accelerated in the channel and their jump- 
ing amplitude increases after each ejection until a max- 
imum is reached as seen in Fig. 2. The more particles 
are injected the smaller is this final amplitude. Beyond 
a certain number n s of particle streams, the trajectories 
however start to lose energy and the overall flux is re- 
duced. This critical value n s characterizes the saturated 
flux q s through Eq. 

In Fig. 3 we see the plot of q s as function of the wind 
velocity it*. Clearly, there exists a critical wind velocity 
threshold it t below which no sand transport occurs at all. 
This agrees well with experimental observations 0, • 
Also shown in Fig. 3 is the best fit to the simulation 
data using the classical expression proposed by Lettau 
and Lettau 4], 
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FIG. 3: Saturated flux q s as function of u*. The full line is 
the fit using the expression proposed by Lettau and Lettau 
Q, 1b oc u*(m* — ut), with ut = 0.35 ± 0.02. The inset shows 
the same data in a double logarithmic plot. The straight 
line corresponds to the least-squares fit to the data of the 
power-law relation, q a oc (u, — u t ) a , with the exponent a — 
2.01 ± 0.01 and the critical point given by u t = 0.33 ± 0.01. 



where C'l is an adjustable parameter. We find very good 
agreement using fit parameters of the same order as those 
of the original work and a threshold value of u t — 0.35 ± 
0.02. This is in fact, to our knowledge, the first time a 
numerical calculation is able to quantitatively reproduce 
this empirical expression and it confirms the validity of 
our simulationprocedure. Other empirical relations from 
the literature [sUEl can a l so be used to fit these results 
[22| . Close to the critical velocity u t interestingly we find 
that a parabolic expression of the form 
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fits the data at least as well as Eq. © , as can be seen in 
the inset of Fig. 3. The threshold obtained in this case, 
namely, u t = 0.33 ± 0.01, is slightly lower than the one 
obtained for the classical cubic expression Eq. ©. We 
believe that our parabolic expression describes correctly 
the critical behavior very close the transition right above 
u t , since at this point the classical assumption of pro- 
portionality between the saltation jump length and u* 
cannot hold. 

Experimentally much more difficult to access is the ve- 
locity profile of the wind within the layer of grain trans- 
port. This profile clearly deviates from the undisturbed 
logarithmic form of Eq. because of the momentum 
the fluid must locally exchange with the particles. In 
Fig. 4 we show the loss of velocity with respect to the 
logarithmic profile without particles of Eq. for differ- 
ent values of q as function of the height y. As clearly 
seen in Fig. 4, the loss of velocity is maximal at the same 
height y m ax, regardless of the value of flux q. Except for 
large values of the flux, dividing the velocity axis by q 



FIG. 4: Profile of the velocity difference u x (0) — u x (q) for 
different values of the flux q at u, = 0.51. The inset shows the 
data collapse of these data obtained by rescaling the velocity 
difference with the corresponding q. 
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FIG. 5: Height y ma x of the maximum loss of velocity as func- 
tion of u*. The height y s of the largest probability to find 
a particle coincides with y m ax- The solid line corresponds to 
the best linear fit to the data with a slope equal to 0.35. By 
extrapolation, the intercept with the 3>axis provides an alter- 
native estimate for the critical point, Ut = 0.35 m/s, that is 
consistent with the other calculations. 



one can collapse all the profiles quite well on top of each 
other as can be verified in the inset of Fig. 4. 

The position y max of the height of maximum loss de- 
pends essentially linearly on as shown in Fig. 5. This 
is consistent with the observation that the saltation jump 
length is proportional to it* . The proportionality con- 
stant obtained from the best linear fit to the data is 
0.35 s. Quantitatively the data in Fig. 5 also fit very well 
into the experiment data plots of Ref. 0] . By extrapo- 
lation to y m ax = 0, we obtain an alternative estimate for 
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the threshold velocity, u t — 0.35 m/s, that is consistent 
with the values calculated before from the fits to the data 
using Eqs. © and Q. 

Whoever has been in the desert or on a beach during a 
very windy day knows that the saltation process in nature 
looks like a sheet of particles floating above ground at a 
certain height y s which strongly depends on the wind 
velocity. This height corresponds to the position of the 
largest likelihood to find a particle as obtained from the 
maximum of the density profile of particles as function 
of height y. Fig. 5 implies that the profile of velocity 
difference of the wind has a minimum at a similar height, 
which is consistent with the maximal loss of momentum. 
Within the error bars our results in fact yield that y s 
coincides with the values of Umax in Fig- 5. More details 
will be presented in Ref. |23. It is important to note 
that both heights, y m ax and y s , also have the same linear 
dependence on u*. 

We have shown in this letter results of simulations giv- 
ing insight about the layer of granular transport under 
conditions of turbulent flow. The lack of experimental 



noise allows for a precise study close to the critical thresh- 
old velocity u t that lead us to a parabolic dependence of 
the saturated flux. In addition, we show that the ve- 
locity profile disturbed by the presence of grains scales 
linearly with the flux of grains, except close to satura- 
tion. Notably a characteristic height appears at which 
the momentum loss in the fluid and the grain density are 
maximized. Moreover, this height increases linearly with 
the wind velocity m* . It would be very interesting to ver- 
ify experimentally these novel predictions. The present 
model can be extended in many ways including the study 
of the dependence of the aeolian transport layer on the 
grain diameter, the gas viscosity, and the solid or fluid 
densities. This would allow to calculate, for instance, 
the granular transport on Mars and compare with the 
expression presented in Ref. [lOj. Work in this direction 
is under way p^ . 
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